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Abstract 

In this paper, the method of discrete exterior calculus for numer- 
ically solving Maxwell's equations in space manifold and the time is 
discussed, which is a kind of lattice gauge theory. The analysis of its 
stable condition and error is also accomplished. This algorithm has 
been implemented on C-|--|- plateform for simulating TE/M waves in 
vacuum. 
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1 Introduction 

Computational electromagnetism is concerned with the numerical study of 
Maxwell's equations. The Yee scheme is known as finite difference time 
domain and is one of the most successful numerical methods, particularly 
in the area of microwave problems [1]. It preserves important structural 
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features of Maxwell's equations [2-6]. Bossavit et al present the Yee-like 
scheme and extend Yee scheme to unstructured grids. This scheme combines 
the best attributes of the finite element method (unstructured grids) and 
Yee scheme (preserving geometric structure) [7,8]. Stern et al [9] generalize 
the Yee scheme to unstructured grids not just in space, but in 4-dimensional 
spacetime by discrete exterior calculus(DEC) [10-20]. This relaxes the need 
to take uniform time steps. 

In this paper, we generalize the Yee scheme to the discrete space manifold 
and the time. The spacetime manifold used here is split as a product of ID 
time and 2D or 3D space manifold. The space manifold can be approximated 
by triangular and tetrahedrons depending on dimension, and the time by 
segments. So the spacetime manifold is approximated by prism lattice, on 
which the discrete Lorentz metric can be defined. 

1. With the technique of discrete exterior calculus, the M value discrete 
connection, curvature and Bianchi identity are defined on prim lattice. 
With discrete variation of an inner product of discrete 1— forms and 
their dual forms, the discrete source equation and continuity equation 
are derived. 

2. Those equations compose the discrete Maxwell's equations in vacuum 
case, which just need the local information of triangulated manifold 
such as length and area. 

The discrete Maxwell's equations here can be re-grouped into two sets of ex- 
phcit iterative schemes for TE and TM waves, respectively. Those schemes 
can directly use acute triangular, rectangular, regular polygon and their com- 
bination, which has been implemented on C-|— I- plateform to simulate the 
electromagnetic waves propagation and interference on manifold. 

2 DEC for Maxwell's equations 

Maxwell's equations can be simply expressed once the language of exterior 
differential forms is used. The electric and magnetic fields are jointly de- 
scribed by a curvature 2— form F in a 4-D spacetime manifold. The Maxwell's 
equations reduce to the Bianchi identity and the source equation 

dF^O d*F^*J (1) 
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where d denotes the exterior differential operator, * denotes the Hodge star 
operator, and 1-form J is called the electric current form satisfying the con- 
tinuity equation 

d*J = 0. 

As the exterior derivative is defined on any manifold, the differential 
form version of the Bianchi identity makes sense for any 3D or 4D spacetime 
manifold, whereas the source equation is defined if the manifold is oriented 
and has a Lorentz metric. Now, we introduce the discrete counterpart of 
those differential geometric objects to derive the numerical computational 
schemes for Maxwell's equations. 

Discrete Lorentz metric 

The spacetime manifold used here is split as a product of ID time and 2D 
or 3D space manifold. The 2D or 3D space manifold can be approximated 
by triangular or tetrahedrons, and the time by segments. The length of 
edge and area of triangular and volume of tetrahedrons gives the discrete 
Riemann metric on space grids. The metric on time grid is the minus of 
length square. The spacetime manifold is approximated by prism lattice, on 
which the discrete Lorentz metric can be defined as the product of discrete 
metric on space and time. 

Discrete exterior calculus 

A discrete differential /c-form. A; e Z, is the evaluation of the differential Re- 
form on all /c-simphces. Dual forms, i.e., forms that we evaluate on the dual 
cell. Suppose each simplex contains its circumcenter. The circumcentric dual 
cell -D((To) of simplex (Tq is 

D{cro) ■= U Int(c((To)c(cTi) • • • c((t^)), 

where (Tj is all the simplices which contains ctq,..., crj_i, and c(a"j) is the 
circumcenter of (Jj. 

The two operators in Eqs.(l) can be discretized as follows: 

1. Discrete exterior differential operator d, this operator is the transpose 
of the incidence matrix of /c-cells on A; + 1-cells. 
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2. Discrete Hodge Star *, the operator scales the cells by the volumes of 
the corresponding dual and primal cells. 

Discrete connection and curvature 

Discrete connection 1— form or gauge field A assigns to each element in the 
set of edges E an element of the gauge group M: 

A: E ^R. 

Discrete curvature 2— form is the discrete exterior derivative of the discrete 
connection 1— form 

F = dA: P^R. 

The value of F on each element in the set of triangular P is the coefficient 
of Holonomy group of this face. The 2— form F automatically satisfies the 
discrete Bianchi identity 

dF = 0. (2) 

Note that since the gauge group R used here is abehan, we need not pick 
a starting vertex for the loop. We may traverse the edges in any order, so 
long as taking orientations into account. 

Discrete Maxwell's equations 

For source case, we need discrete current 1— form J. Let A 
Lagrangian functional be 

L{A,J) = -l{dA,dA) + {A,J), 

where 

{dA,dA) := (^)lx|E|(rf)|E|x|F|(*)|F|x|F|(d)^|x|E|(^)|^|xl 
{A, J) := {A)i,^\E\{*)\E\x\E\{j'^)\E\xl 

Supposing that there is a variation of Ai, vanishing on the boundary, we 



— ^Ai and the 

E 
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have 

OaMAJ) = dAA-W^^dA) + {A,J)) 

= dAi{-^{A)iy,iE\{d)\E\x\F\i*)\F\x\F\{d)fp\^\E\{A)fj^^^-^ 

+ i^)lx\E\i*)\E\x\E\{J^)\E\xl) 

= ■•■'^J^' ^)lx\E\{d)\E\x\F\{*)\F\x\F\{d)fF\^\E\i^)fE\xl 

i 

-|(^)lx|i?|(c^)|E|x|F|(*)|F|x|F|((^)j^|x|E|(0,.-,^l^, •••,0)j^,|xl 

i 

+ (0) •••'s^^' •••'0)lx|B|(*)|E|x|£;|('^'^)|E|xl 

j 

= -(0> ■••'0)lx|E|(c?)|S|x|F|(*)|F|x|F|(c?)jS'|x|S|(^)js|xl 

|E|(*)|E|x|E|(-^^)|E|) 



+ (0, •••,0)ix|E|(*)|E|x|E|(-'^^)|E|xl 



The Hamilton's principle of stationary action states that this variation 
must equal zero for any such vary of Ai, implying the Euler-Lagrange equa- 
tions 

-(0?)|E|x|F|(*)|F|x|F|(0?)jS'|x|E|(^)|E|xl + (*)|i=;|x|E|(./'^)|E|xl = 0, 

which is the discrete source equation 

5F = J, (3) 

where S = Since {d^y = 0, the discrete continuity equation can 

express as: 

(f*J=0. (4) 

The equations of discrete Bianchi identity (2), source equation (3), and con- 
tinuity equation (4) are called discrete Maxwell's equations. 

Discrete Gauge transformations 

Discrete gauge transformations are maps 

A^A + df 

for any 0— form or scalar function / on vertex. Since the discrete exterior 
derivative maps 

F ^ F + d^f ^ F, 
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the discrete Maxwell's equations(2-4) are invariant under discrete gauge trans- 
formations. Since the discrete continuity equation (4) ensures 



L{A + df,J) = -'^{d{A + df),d{A + df)) + {A + df,J)^L{A,J). 



That is to say the Lagrangian function is also invariant under discrete gauge 
transformations . 

3 Explicit schemes 
Schemes for TE wave 

The discrete current 1— form, discrete curvature 2— form and its dual can be 
written as 

J = i-pedt, Je) = A dt + B"" * = i/" A rft - D"-^ , 

where n and n + ^ denote the coordinate of the time, E = EiC^ (electric 

E 

field) is discrete 1— form on space, B — ^ BiP^ (magnetic field) is discrete 



2— form on space, H — Y^Hi* P'^ (magnetizing field) is the dual of B on 



space, D — Di * (electric displacement field) is the dual of E on space, 

E 

Pedt (charge density) is the discrete 1— form on time, Je = XI Jei^^ (electric 

E 

current density) is the discrete 1— form on space. The discrete Maxwell's 
equations can be rewritten as 



{df, J) = (/)ix|y|(cF)|y|x|£;|(*)|£;|x|E|(J^)|E|xi = 0, 



so we have 



p 



p 



dsB^ 

dsE'^+h Adt = 
djD^--2 
d^H" Adt 





-dtB"" 

*{Pedtr-'2 



where ds, are the restriction of d, d^ 



on space, and 



dtB 



Qn+l _ 



A dt djD^-^ : 



Adt. 



(5) 



in 



At 



At 
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If allowing for the possibility of magnetic charges and current discrete 3— form 

J = {Pm, -Jm A dt), 
the symmetric scheme can be written as 



(6) 



dsB^ = p^^ 

dsE^+^Adt = -dtB"" - J'!^'^ Adt 
d^^D'^-l = *{pedtY-^ 
djH^'Adt = djD''-^ + *J2, 

where 

Pm — ^ Pmi^' (ma-giietic charges) is discrete 3-form on space, 

Tet 

Jm — Yli JmiP^ (current) is discrete 2— form on space. 

p 

The compact form of Eqs. (6) can be written as 

with discrete continuity equations or integrability conditions 

dJ = (f*J=0. 

Proposition 3.1 // the initial condition satisfies the first and third equa- 
tions in Eqs. (6), the solution of the second and fourth equations in Eqs. (6) 
automatically satisfy Eqs. (6). 

Proof. Because the dimension of spacetime is 3 + 1, therefore 

d^ * (pedt) —0 dj * Je — dsPm — dtJm A dt — 0., 
and the continuity equations can be reduced to 

dJi^pedt""-^) -d'^*J2^0 - dsf^'^ Adt + dtpl = 0. 
So we have 

rff -d[* {pedtY-^2 ^ -d[* {pJtY''^ - Adt- *J^) 

= 

dtdsB^-dtpl = -dtp';^ + d,{dsE^+-2 Adt + jZ^~' Adt) 

= 0. 
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□ 

Now we show the scheme (5) on the product of 2D discrete space manifold 
and time. The second and fourth equations in Eqs.(5) based on Fig.l are 



n-i 

2 n 2 



At 



I jn 

1 



* Ci 





ei 




62 






At 


Pi 





(7) 



where | | denotes the measure of forms and dual. The summation on the 
right is orient, that is to say, inverse the orientation of Cj, then multiply — 1 
with Ei. Eqs.(7) can be implemented on 2D discrete manifold directly(see 
Fig.l). Eq.(7) on rectangular gird is just the Yee scheme. 




Figure 1: edge and face with direction 
In the absence of magnetic or dielectric materials, the relations are simple: 

A = ^oEi Bi = fioHi, (8) 

where Eq and /io are two universal constants, called the permittivity of free 
space and permeability of free space, respectively. With relations (8), Eqs.(7) 
can be rewritten into an explicit iterative scheme for TE wave. 

jj^r^ - Er'' _^ jn _ - 

^0 — <~ ^el — 
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At 
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(9) 
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The symmetric TE wave scheme induced from Eqs.(6) can be written as 
follows. 



At 



el 



'At 





1 1 1 

ei + -C/2 


62! + £^3 ^lesl 


1^1 





> TE (10) 



Schemes for TM wave 

If writing 



= i/^+l Adt-D''' 

J = {-pe,Je/\dt) 



= -S" Adt- S"-2 

J ( Pmdt, Jrn)i 



where H — '^HiC^ is the discrete 1— form on space, D — J^^iP^ is the 



E P 

discrete 2— form on space, E = J2 * is the dual of D on space, B = 

*p 

Y^Bi* is the dual of H on space, pe = PeiT^ is the discrete 3— form 

*E Tet 

on space, Jg = Yl JeiP^ is the discrete 2— form on space, Pmdt is the discrete 

p 

1— form on time, Jm — Yl JmiE"^ is the discrete 1— form on space, the discrete 

E 

Maxwell's equations can be rewritten as 



P'e 



dsH'''+^ Adt ^ dtD"" + Je^^ A dt 



diB 
dTE"" A dt 



= *{pmdtY 2 



(11) 



Proposition 3.2 // the initial condition satisfies the first and third equa- 
tions in Eqs.(ll), the solution of the second and fourth equations in Eqs.(ll) 
automatically satisfy Eqs.(ll). 

Proof. Because the dimension of spacetime is 3 + 1 or 2 + 1, therefore 
* (Prndt) — dj * Jm — dgPe — dfJe A dt — 0, 
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and the continuity equations can be reduced to 

dj * {pmdty-^ - dl{*Jl) = 4 Je"^^ ^dt- dtpl = 0. 

So we have 
dtdsD^ - dtp'l 



= -dtpl - dsidsH^'+l Adt- Je^^ A dt) 
= 

d'^d^B''-^ -df* {pmdtY-^ = -df * {pmdtY-^ + d:^{d^E^+^ Adt + *J^) 

= 0. 

□ 

Now we show the scheme (11) on the product of 2D discrete space mani- 
fold and time. The second and fourth equations in Eqs.(ll) based on Fig.l 
are 



At 



■\n+l r-tn 



* Cl 



At 





ei 


+ <^ 


62 
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Pi 









(12) 



With relations (8), Eqs.(12) can be rewritten into an explicit iterative scheme 
for TM wave. 





ei 




62 













A^o —J. r Jjni - - 



> TM (13) 



At 



General schemes 

For real world materials, the constitutive relations are not simple propor- 
tionalities, except approximately. The relations can usually still be written: 

D^eE iiH, 

but e and p are not, in general, simple constants, but rather functions. With 
Ohm's law 



E — — J, J„ 



-H, 
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where a is the electrical conductivity and cr^ is magnetic conductivity. The 
DEC schemes can be written as 



-,n+i ^n- 



E';'^^ - E'; ^ E"^^^ + E'l ^ HI' - 
e h (7 



A/ 2 

n+l un Ui^+i _L 
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At 2 

ffi"^^ - H^~^ H^^^ + i/""^ E^ - E^ 

^ X7 ^ ^"^ n = \ ]~ 

At 2 * ^1 



> TM 



4 Stability, convergence and accuracy 
Stability 

The Courant-Friedrichs-Lewy condition is a necessary condition for conver- 
gence while solving certain partial differential equations numerically. Now, 
we find this condition for scheme (10). Condition for scheme (13) can be 
induced in the same way. First, we decompose DEC algorithm into temporal 
and spacial eigenvalue problems. 
The temporal eigenvalue problem: 



It can approximated by difference equation 



{Aty 



A^o- (14) 



Supposing 



and 



= h;^ cos(niAt) H^-^ = cos(n2At) 
H^+^ = sin(ni At) H^'^ = sin(n2At), 
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and substituting those into Eq.(14), we obtain 

cos(r;,iAt) + cos(n2At) — 2 

sin(?T,iAt) + sin(n2At) — 2 



therefore 



(At)2 
4 



A, 



A, 



< A < 0. 



This is the stabile condition for the temporal eigenvalue problem. 
The spacial eigenvalue problem: 

c^AH = AH 

It can be approximated by difference equation (15) based on Fig. 2. 

C tAO 'BO 'CO 




Let 



Figure 2: Face and dual face 



Hi = Ho cos{doi) or H, = Hq sm{doi) , 
and substitute into Eq.(15) to obtain 

^^A = ^(cos(c/oyi) - 1) + ^(cos(c/oi3) " 1) + 7^(cos(c/oc) 



hi , 



CO 



^^A = ^(sin(c/oA) - 1) + ^(sin(doB) - 1) + ^(sin(doc; 

C t^O l-BO 'CO 
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So we have 



2c2 



Pi 



123 



^ + ^ + ^1 <A<0. 

I- AO I'BO ho. 



In order to keep the stability of scheme (13), we need 



or 



2< 



(At) 



At < Mimp^23ep 



-Pl23 
/ 



^23 _|_ ^34 _|_ ^45 

Iao I'BO Ico 



2P 



123 



^23 _|_ ^34 ^45 
\ V^O ^BO ho J ) 



(16) 



Convergence 

By the definition of truncation error, the exact solution of Maxwell's equa- 
tions satisfy the same relation as DEC scheme except for an additional term 
0{{At)^ + At| * e|). This expresses the consistency, and so convergence for 
DEC scheme by Lax equivalence theorem (consistency + stability = conver- 
gence) . 



Accuracy 

The derivative of Maxwell's equations is approximated by first order differ- 
ence in schemes (10) and (13). Equivalently, H and E are approximated 
by linear interpolation functions. Consulting the definition about accuracy 
of finite volume method, we can also say that schemes (10) and (13) have 
first order spacial and temporal accuracy, and have second order spacial and 
temporal accuracy on rectangular grid with equivalent space and time steps. 
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5 Implementation 

The DEC algorithm of Maxwell's equations was implemented in C++ plat- 
form. The Fig. 3 shows the flowchart of DEC schemes for Maxwell's equations. 
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Figure 3 

In the common practice, not every simulation step needs to be visualized, 
especially when the time step size is too small. Fig. 4 exhibit Gaussian pluses' 
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waveforms simulated by DEC. 




Figure 4: Simulation of Gaussian pluse on Stanford bunny by DEC 

Fig. 5 exhibits two sources Gaussian pluses interference simulated by DEC. 
Our algorithm can simulate more complex situation on surface and 3D-space 
manifold. 




Figure 5: Simulation of Gaussian pluse interference on sphere by DEC 
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